% fig3 分岔图和李指数图(Bifurcation&LEs)
clear; clc; close all;

%% —— 公共参数 —— 
g11 = 1; g12 = 1; g21 = -2; g22 = 2;   % 连接权重
Niter  = 500;   % 总迭代步数
Ntrans = 400;   % 舍弃热身步数

%% —— (a) μ–y 分岔图 设置 —— 
mu_a   = linspace(0,1,600);           % μ 采样点
init_a = [ 0.1, -0.1;                % 两组初始值
          -0.1,  0.1 ];
col_a  = [0.5,1,0.5;                  % 草绿色
         0,1,1];                     % 青色
nMu_a = numel(mu_a);
nIV   = size(init_a,1);

% 预分配并行结果
y_collect = cell(nIV, nMu_a);

parfor i = 1:nMu_a
    mu = mu_a(i);
    tmp = cell(nIV,1);
    for j = 1:nIV
        x = init_a(j,1);
        y = init_a(j,2);
        ys = zeros(Niter - Ntrans,1);
        idx = 0;
        for k = 1:Niter
            x_next = mu*x + g11*sin(x) + g12*sin(y);
            y_next = mu*y + g21*sin(x) + g22*sin(y);
            x = x_next;  y = y_next;
            if k > Ntrans
                idx = idx + 1;
                ys(idx) = y;
            end
        end
        tmp{j} = ys;
    end
    y_collect(:,i) = tmp;
end

%% —— (b) μ–LE 光谱 设置 —— 
mu_b = mu_a;               % 与分岔图同分辨率
x0 = 0.1;  y0 = -0.1;       % LE 初始条件
nMu_b = numel(mu_b);
LE1 = zeros(nMu_b,1);
LE2 = zeros(nMu_b,1);

parfor idx = 1:nMu_b
    mu = mu_b(idx);
    x = x0; y = y0;
    v1 = [1;0]; v2 = [0;1];
    s1 = 0;     s2 = 0;
    for k = 1:Niter
        x_next = mu*x + g11*sin(x) + g12*sin(y);
        y_next = mu*y + g21*sin(x) + g22*sin(y);
        J = [mu + g11*cos(x),  g12*cos(y);
             g21*cos(x),      mu + g22*cos(y)];
        % Gram–Schmidt
        w1 = J*v1; norm1 = norm(w1); v1 = w1/norm1;
        w2 = J*v2 - (v1'*(J*v2))*v1; norm2 = norm(w2); v2 = w2/norm2;
        if k > Ntrans
            s1 = s1 + log(norm1);
            s2 = s2 + log(norm2);
        end
        x = x_next;  y = y_next;
    end
    M = Niter - Ntrans;
    LE1(idx) = s1 / M;
    LE2(idx) = s2 / M;
end

%% —— 绘图 —— 
figure;

% (a) μ–y 分岔图
subplot(1,2,1); hold on;
for i = 1:nMu_a
    mu = mu_a(i);
    for j = 1:nIV
        plot(mu*ones(size(y_collect{j,i})), y_collect{j,i}, ...
             '.', 'Color', col_a(j,:), 'MarkerSize', 1);
    end
end
xlabel('\mu','FontSize',14);
ylabel('y','FontSize',14);
title('(a) \mu–y 分岔图','FontSize',16);
xlim([0,1]); grid on;
legend('(0.1,-0.1)','(-0.1,0.1)','Location','northeast');

% (b) μ–LE 光谱
subplot(1,2,2); hold on;
plot(mu_b, LE1, 'r-', 'LineWidth',1.2);
plot(mu_b, LE2, 'b-', 'LineWidth',1.2);
xlabel('\mu','FontSize',14);
ylabel('LEs','FontSize',14);
title('(b) LE_1 & LE_2 光谱','FontSize',16);
xlim([0,1]); grid on;
legend('LE_1','LE_2','Location','northeast');

hold off;


